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ABSTRACT 


The  Australian  Challenge  was  a  set  of  experiments  chosen  as  a  test  case  to 
establish  how  well  computer  calculations  can  model  the  blast  overpressure 
around  a  complex  set  of  structures.  This  paper  presents  the  first  three 
dimensional  calculation  on  the  Australian  Challenge  using  CFD  codes.  These 
codes  have  been  developed  at  AMRL  and  the  Australian  Challenge  represents 
their  first  test  in  a  'real  world'  scenario.  The  calculations  show  good  agreement 
with  the  experimental  results.  The  program  can  now  be  used  with  a  good  degree 
of  confidence  in  modelling  other  complex  scenarios,  reducing  the  need  for 
expensive  experimentation.  The  safety  aspects  of  any  changes  in  operational 
procedure  can  be  evaluated  quickly  and  efficiently. 
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EXECUTIVE  SUMMARY 


The  1994  Australian  white  paper  on  Defence  "Defending  Australia"  requires  that 
the  Australian  Defence  Forces  be  ready  to  defend  against  a  short  warning 
conflict.  A  short  warning  conflict  could  vary  between  a  series  of  low  level  raids 
to  a  larger,  more  protracted  operation.  In  these  scenarios  it  is  very  conceivable 
that  the  forces  would  require  rapid  entry  into  structures,  usually  buildings,  to 
achieve  their  mission.  The  obvious  way  to  do  this  is  through  an  explosive 
charge  either  on  a  wall,  door  or  window.  The  question  then  becomes  what  is  the 
distance  to  stand  away  from  the  explosive  charge  to  prevent  any  injuries  and 
also  what  is  the  likely  effect  on  the  people  inside  the  building.  This  can  be 
answered  by  performing  an  expensive  series  of  experimental  trials  or  by  using  a 
reliable  computer  simulation.  This  paper  presents  numerical  simulation  for  a 
specific  set  of  scenarios  for  which  experimental  results  are  available  as  a  way  of 
validating  the  computer  program  in  a  realistic  complex  situation. 

The  scenarios  modelled  were  issued  as  a  challenge  to  the  numerical  modellers 
within  the  TTCP  nations.  This  paper  presents  the  first  results  for  these  scenarios 
based  on  computational  fluid  dynamics  (CFD)  using  full  three-dimensional 
geometry.  The  scenarios  consist  of  explosive  charges  detonating  on  a  concrete 
reinforc^  wall  of  a  building  which  may  or  may  not  have  open  windows  and 
doors.  The  blast  overpressure  from  the  explosion  was  calculated  and  compared 
against  the  experimental  results  in  considerable  detail.  This  allowed  for  an 
analysis  of  the  approximations  used  within  the  calculations  and  an  estimate  of 
their  impact  on  the  computed  results. 

The  computational  results  show  that  the  program  has  the  ability  to 
successfully  calculate  the  blast  overpressure  from  an  explosion  in  a  structurally 
complex  scenario.  This  validation  of  the  program  will  allow  it  to  be  used  with 
confidence  to  model  complex  operational  scenarios  requiring  the  rapid  entry  of 
forces  into  buildings.  It  will  reduce  dramatically  the  need  for  expensive  and 
time-consuming  trials  in  this  area  and  allow  hypothetical  situations  to  be 
explored.  It  will  also  allow  the  safety  of  current  operational  procedures  to  be 
evaluated,  as  well  as  any  proposed  changes  to  these  procedures. 


Author 


G.  Kemister 

Weapons  Systems  Division 


Gary  Kemister  graduated  from  Sydney  University  in 
1980  with  a  B.Sc  (Hons  I)  and  in  1985  with  a  Ph.D.  in 
theoretical  chemistry.  After  working  at  Oxford 
University,  University  of  North  Carolina  at  Chapel 
Hill,  Sydney  University  and  La  Trobe  University  on 
various  aspects  of  electronic  structure  of  materials  he 
joined  the  Defence  Department  in  1991  working  for 
the  Analytical  Studies  Unit  in  Force  Development  and 
Analysis  Division  in  Canberra.  In  1993  he  moved  to 
Explosive  Ordnance  Division  (now  part  of  Weapons 
Systems  Division)  at  Maribyrnong  and  is  currently 
working  on  numerical  simulations  of  the  physics  and 
chemistry  of  explosions. 


Contents 


1.  INTRODUCTION . 1 

2.  EXPERIMENTAL  DETAILS . 1 

3.  COMPUTATIONAL  DETAILS . 6 

4.  RESULTS . 8 

5.  CONCLUSIONS . 23 

6.  ACKNOWLEDGEMENTS . 24 

7.  REFERENCES . 25 


APPENDIX  1:  Inclusion  of  an  Arbitrary  Geometry  in  the  Program 


27 


DSTO-TR-0224 


1.  Introduction 


The  1994  Australian  white  paper  on  Defence  "Defending  Australia"  requires  that  the 
Australian  Defence  Forces  be  ready  to  defend  against  a  short  warning  conflict.  A  short 
warning  conflict  could  vary  between  a  series  of  low  level  raids  to  a  larger,  more 
protracted  operation.  In  these  scenarios  it  is  very  conceivable  that  the  forces  would 
require  rapid  entry  into  structures,  usually  buildings,  to  achieve  their  mission.  The 
obvious  way  to  do  this  is  through  an  explosive  charge  either  on  a  wall,  door  or 
window.  The  question  then  becomes  what  is  the  distance  to  stand  away  from  the 
explosive  charge  to  prevent  any  injuries  and  also  what  is  the  likely  effect  on  the  people 
inside  the  building.  To  answer  these  questions  experimental  trials  have  been 
conducted  in  selected  scenarios  and  the  overpressures  measured.  However 
experimental  trials  are  expensive  and  specific  to  the  selected  scenario,  consequently  it 
is  better  to  be  able  to  model  scenarios  and  their  variations  with  a  computer  simulation 
which  has  been  validated  against  experimental  measurements.  This  work  is  the 
validation  of  the  computer  model  developed  at  WSD  AMRL  against  a  set  of 
experimental  overpressure  measurements  in  a  complex  building  scenario. 

The  Australian  Challenge  consists  of  two  sets  of  air  blast  data  recorded  by  AMRL 
in  1990  at  Woomera,  South  Australia,  as  part  of  the  collaborative  demolition  trials 
conducted  under  the  auspices  of  The  Technical  Cooperation  Program  (TTCP)  [1].  The 
test  series  from  which  the  results  were  generated  comprised  a  number  of  firings  of 
experimental  contact  breaching  charges  against  the  walls  of  derelict  reinforced 
concrete  buildings.  Air  blast  gauges  were  positioned  around  the  building  returns  to 
assess  the  feasibility  of  personnel  firing  the  charges  in  close  proximity,  but  out  of  line 
of  sight.  The  air  blast  results  have  been  circulated  around  the  TTCP  nations  as  a  test 
for  the  modelling  capability  of  blast  overpressure  calculations.  Several  attempts  [2,3,4] 
have  been  made  to  model  these  experimental  results  with  generally  good  agreement, 
however  the  previous  calculations  have  been  limited  by  approximations  concerning 
the  geometry  of  the  experimental  set-up.  In  the  calculations  reported  here  we  have 
modelled  the  buildings  in  full  3-D  representation  and  compared  the  results  to  both  the 
experimental  and  previous  computational  results. 


2.  Experimental  Details 


The  Australian  Challenge  consists  of  sets  of  experimental  results  for  two  selected 
building  layouts.  The  layout  for  the  first  experiment,  event  number  three,  is  shown  in 
Figure  1  and  the  layout  for  the  second  experiment,  event  number  ten,  is  shown  in 
Figure  2.  Note  that  event  3  and  event  10  were  conducted  at  opposite  ends  of  the 
building  structures.  Photographs  of  the  buildings  showing  damage  sustained  prior  to 
the  conduct  of  the  tests  being  analysed  in  this  paper  are  shown  in  Figure  3.  These 
photographs  are  taken  from  the  end  of  the  buildings  where  event  3  was  conducted. 


SHOT  NUMBER!  3 
DATD12  JULY  90 
TIME'  10*10 
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Figure  2.  Experimental  layout  for  event  number  ten. 
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In  each  experiment  a  charge  was  detonated  to  breach  the  wall  of  the  building  and 
blast  overpressure  measurements  were  recorded  at  four  locations  outside  the  building. 
The  breaching  charges,  in  each  case,  consisted  of  four  copper  lined  linear  shaped 
charges,  each  500mm  in  length,  arranged  in  a  square  array.  Plastic  Explosive  No.  4 
(88%  RDX,  11%  grease  and  1%  PEDO)  was  used  to  fill  the  charges  which  were 
initiated  simultaneously  with  detonating  cord  (Fig  4). 


Figure  4.  Sketch  of  demolition  charge. 

Event  3  -  For  event  3  the  charge  weight  was  8.72  kg.  Event  3  is  characterised  by  a  free 
standing  wall  (1.8m  high)  with  a  bend  in  it  around  the  comer  from  where  the  charge 
is  detonated  (see  Figure  1). 
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Event  10  -  For  event  10  the  charge  weight  was  21.6  kg.  In  event  10  the  building  is 
characterised  by  two  holes  caused  by  previous  events  which  for  our  purposes  will  be 
considered  as  windows  and  a  door.  There  is  also  a  free  standing  wall  (1.8m  high) 
parallel  to  the  structure  (see  Figure  2). 

Pressure  gauges  were  placed  at  the  positions  PI,  P2,  P3  and  P4  marked  on  Figures 
1  and  2.  The  experimental  time-pressure  traces  were  recorded  for  each  gauge  and 
each  event  at  several  sampling  speeds,  for  example,  the  time-pressure  trace  for  gauge 
1  in  event  3  was  recorded  at  time  intervals  of  8  ns,  40  ns  and  160  ns  with  each  being 
recorded  for  512  steps.  Thus  the  recordings  span  different  time  lengths,  being  4.096 
ms,  20.48  ms  and  81.92  ms  respectively.  All  of  the  experimental  traces  used  in  this 
paper  are  for  the  longest  time  step  i.e.,  160  ns,  extending  out  to  ~80  ms.  These  traces 
were  chosen  because  they  most  closely  correspond  to  the  time  step  used  in  the 
calculations.  The  time  step  is  initially  80  ns  but  increases  in  length  as  the  speed  of  the 
pressure  front  decreases.  Variable  time  steps  were  used  to  improve  the  efficiency  of 
the  calculation  as  they  can  cut  the  number  of  time  steps  requir^  to  reach  a  particular 
time  point  by  approximately  one-third.  For  example  the  time  step  is  about  320  ns  by 
the  time  the  calculation  has  finished. 

Further  experimental  details  can  be  obtained  from  Reference  1. 


3.  Computational  Details 


The  calculations  were  performed  using  a  3-dimensional  finite  difference  computer 
program,  developed  at  AMRL  [5],  which  solves  the  Euler  equations  for  the 
conservation  of  mass,  momentum  and  energy  of  an  inviscid,  compressible  fluid  (in 
this  case  air).  A  3-dimensional  Cartesian  grid  with  fixed  spacing  was  used  as  a  basis 
for  the  calculation.  Operator  splitting,  which  reduces  the  complexity  of  a  3- 
dimensional  calculation  down  to  three  1-dimensional  calculations  at  each  time-step 
and  grid  point,  was  also  used.  A  1-dimensional  Flux-Corrected  Transport  (FCT) 
algorithm,  with  fourth  order  phase  accxiracy  and  an  overall  second  order  accuracy  on 
uniform  grids,  was  used  to  solve  the  1-dimensional  calculations  [6].  Rigid  wall 
boundary  conditions  were  used  for  the  walls  and  the  ground  in  the  calculations  whilst 
outflow  boundary  conditions  were  used  at  all  the  other  edges  of  the  computational 
grid.  The  computational  cell  size  was  5x  =  6y  =  5z  =  20  cm  and  the  largest  grid  size 
used  was  95  x  130  x  60  giving  741,000  cells.  The  time  step  was  calculated  at  each 
iteration.  In  order  to  perform  the  calculations  without  reprogramming  the  code  each 
time  to  take  into  accoimt  different  geometries,  modifications  were  made  to  the  code  to 
allow  the  input  of  an  arbitrary  shape  (within  the  restriction  of  the  fixed  C!artesian  grid) 
and  these  modifications  are  described  in  more  detail  in  Appendix  1.  The  calculations 
were  performed  on  either  the  EOD(S)  RISC  computer  (eodrisc-6000s.dsto.gov.au)  or 
the  EOD(M)  HP  workstation  (eod-hpl. dsto.gov.au). 
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To  begin  the  caloilations  requires  the  simulation  of  an  explosion  on  the  wall.  This 
can  be  modelled  using  the  full  chemistry  and  physics  of  the  explosion;  however  this  is 
both  difficult  and  time-consuming  in  a  calculation  of  this  size.  Furthermore  it  contains 
a  large  amount  of  detailed  information  close  in  to  the  explosion  that  is  simply  washed 
out  by  the  time  the  pressure  pulses  reach  any  of  the  pressure  gauges.  Combustion 
Dynamics  [3]  have  developed  an  approximate  method  of  starting  calculations  without 
modelling  the  full  explosion  which  gives  the  correct  behaviour  in  the  far-field,  i.e.,  a 
sufficient  distance  away  from  the  explosion.  In  this  approximation  1kg  of  TNT  (or 
equivalent)  is  replaced  by  0.0164  m^  of  gas  at  a  pressure  of  1000  atmospheres  and  a 
temperature  of  3000K.  Thus,  for  example  in  event  3,  8.72  kg  of  PE4  with  a  TNT 
equivalence  of  1.35  gives  an  equivalent  weight  of  TNT  of  11.772  kg  and  a  volume  of 
0.193  m^.  The  gas  bubble  approximation  has  been  used  in  all  the  calculations  reported 
here. 

The  gas  bubble  approximation  to  an  explosive  is  designed  to  give  the  correct 
energy  content  for  the  explosive.  In  the  far-field  this  will  give  the  correct  time- 
pressure  trace  to  within  the  computational  accuracy.  However,  in  the  near  field  it  will 
be  unable  to  accurately  represent  the  sharp  peak  caused  by  the  explosive  because  of 
the  assumption  of  an  equal  pressure  everywhere  inside  the  gas  bubble  and  because  of 
the  finite  size  of  the  grid.  The  finite  grid  size  can  be  overcome  by  making  it 
sufficiently  small  however  this  can  result  in  the  calculation  either  being  too  large  to  fit 
in  the  memory  of  the  available  workstations  or  requiring  an  inordinate  length  of  time 
to  perform  the  calculations.  Another  factor  affecting  accuracy  in  the  near  field  is  the 
shape  of  the  explosive.  The  explosive  has  a  much  smaller  volume  than  the  gas  bubble 
and  therefore  any  fine  details  of  the  explosive  shape  may  be  lost  when  it  is  replaced  by 
the  larger  gas  bubble  volume  which  in  itself  is  restricted  to  cells  defined  by  the  grid. 
Whether  these  details  are  important  depends  on  the  distance  at  which  the  interactions 
with  the  structures  occur  and  at  which  the  measurements  are  taken. 

The  four  linear  shaped  charges  used  in  the  experiments  were  arranged  in  a  square 
array  as  shown  in  Fig  4.  Given  the  grid  size  that  is  needed  to  make  the  full 
computation  feasible  within  a  reasonable  length  of  time  (i.e.,  20  cm  cell  size),  the 
explosive  is  approximated  by  a  volume  of  gas  1  cell  deep  (20cm)  by  2  cells  wide 
(40  cm)  by  the  length  required  to  make  the  correct  volume  of  the  equivalent  explosive 
weight,  diagrammatically  shown  in  Figure  5.  The  size  of  the  square  is  enlarged  to 
avoid  any  problems  caused  by  overlapping  squares  and  to  maintain  the  correct 
relative  geometry.  This  is  a  coarse  approximation  in  close  but  it  is  reasonable  in  the 
far-field.  One  way  of  improving  this  approximation  is  to  break  the  calculation  up  into 
separate  problems,  that  is,  near  field  measurements  and  far-field  measurements, 
where  a  much  finer  grid  is  used  to  calculate  the  near  field  measurements  (e.g.,  2  cm). 
This  will  allow  more  cells  in  which  to  resolve  the  peak  pressure  in  the  near  field  at  the 
expense  of  doing  an  extra  calculation. 
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Figure  5.  Diagrammatic  sketch  of  the  explosive  setup.  The  solid  lines  are  the  explosive  and  the 
dashed  lines  are  the  'gas  bubble'  approximation. 


4.  Results 


Combustion  Dynamics  have  modelled  the  Aiistralian  Challenge  in  a  two  dimensional 
axisymmetric  grid  [3].  In  this  approximation  the  axis  of  symmetry  is  through  the 
centre  of  the  charge  on  the  wall  and  everything  is  then  radially  symmetric  about  that 
axis.  Thus  the  free  standing  wall  is  modelled  as  a  cylinder  with  a  radius  of 
approximately  10m  from  the  charge.  Their  results  show  good  agreement  for  the  first 
peak  reaching  the  pressure  gauge  but  the  subsequent  peaks  are  modelled  less 
accurately  due  to  the  inability  to  model  the  1.8  m  wall  correctly  in  the  cylindrical 
coordinates. 

Jones  [2]  has  modelled  event  3  of  the  Australian  Challenge  in  a  similar  way  to  this 
work  but  only  looking  at  the  first  3  gauges  due  to  computational  constraints.  The 
program  used  in  this  paper  is  a  development  of  the  Jones  program. 

There  was  no  time  of  arrival  data  available  from  the  experiments  and  so,  for 
comparison  purposes,  the  initial  arrival  of  the  experimental  curve  and  the 
computational  curve  have  been  superimposed.  This  allows  for  analysis  of  the  relative 
time  of  arrival  of  the  various  peaks. 

Event  3. 


The  experimental  results  (solid  line)  and  the  results  of  our  calculations  (dashed  line) 
for  pressure  gauges  1  to  4  in  Event  3  are  given  in  Figures  6,  8,  9,  and  12  respectively. 
The  peak  pressures  for  this  event  are  tabulated  in  Table  1. 
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Table  1 .  Peak  experimental  and  calculational  pressures  for  gauges  1  to  4  in  Event  3. 


Gauge  Number 

Peak  Overpressure  (kPa) 

Experiment 

Peak  Overpressure  (kPa) 

Simulation 

1 

1130 

636 

2 

47 

42 

3 

11.7 

15 

4 

5.0 

10 

It  can  be  seen  from  Figure  6  and  especially  Table  1  that  the  peak  pressure  given  by 
the  computations  at  gauge  1  is  significantly  lower  than  the  experimentally  recorded 
peak  pressure.  The  single  pressure  pulse  with  little  other  structure  seen  in  Figure  6  is 
the  expected  result  for  a  pressure  gauge  at  this  distance  from  this  size  explosion.  The 
small  pulse  in  the  calculation  at  ~  20  ms  is  the  secondary  expansion  of  the  explosive. 
The  experimental  pulse  for  the  secondary  expansion  is  not  obvious  in  Figure  6, 
presumably  being  swamped  by  the  large  effects  of  the  main  peak. 


Figure  6.  Plot  of  Overpressure  versus  Time  for  Gauge  1  in  Event  3. 
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The  lower  peak  pressure  given  by  the  calculation  is  due  to  the  coarse  grid  size 
(20  cm  cells)  used  in  the  full  calculation.  By  repeating  the  calculation  with  a  grid  size 
of  4  cm  and  a  correspondingly  smaller  timestep,  but  only  looking  at  gauge  1,  (it  is 
impractical  due  to  the  size  of  the  grid  to  include  any  other  gauge  positions),  the  effect 
of  the  ceU  size  on  the  calculation  can  be  inferred.  The  pressure  plot  for  this  calculation 
is  given  in  Figure  7. 


Figure  7.  Plot  of  Overpressure  vs  Time  for  Gauge  1  in  Event  3  using  a  finer  computational 
grid  size  and  timestep  and  comparing  it  against  a  finer  time  recording  of  the  experimental 
results. 

The  experimental  peak  pressure  in  Figure  7  is  1550  kPa  and  is  higher  than  the 
value  quoted  in  Table  1  (1130  kPa)  since  the  experimental  sampling  rate  is  higher  in 
this  reading.  That  is,  for  this  comparison,  the  calculations  are  being  compared  with 
the  experimental  results  which  were  taken  at  8ns  time  intervals  rather  than  the  160  ns 
used  for  the  comparisons  in  the  rest  of  this  paper.  For  this  section  of  the  work  the 
calcvilations  were  performed  with  a  timestep  in  the  range  4-32  ns,  varying  with  the 
speed  of  the  pressure  front  (as  for  the  other  calculations  in  this  paper).  The  peak 
pressure  for  the  calcvilation  in  Figure  7  is  1100  kPa  (compared  to  the  value  of  636  kPa 
given  in  Table  1). 

Figure  7  also  shows  that  the  calculation  does  not  resolve  the  sharpness  of  the  initial 
rise  in  pressure  to  match  the  experiment.  This  is  directly  related  to  the  cell  size  and 
timestep  in  the  calculation.  The  calculation  can  not  hope  to  resolve  a  rise  time  of  the 
order  of  1  ns  if  the  timestep  is  of  the  order  of  4  ns,  thus  the  pressure  curve  from  the 
calculation  in  Figure  7  has  a  much  more  rounded  shape  than  the  experimental  curve. 
Apart  from  the  resolution  of  the  initial  peak  there  is  good  agreement  between  the 
experimental  and  computational  curves. 
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The  pressure  trace  given  in  Figure  8  and  the  peak  pressures  given  in  Table  1  for 
gauge  2  show  quite  good  agreement  between  the  experiment  and  the  calculation.  As 
the  gauge  is  just  around  the  comer  there  is  no  discernible  reflection  from  the  adjacent 
concrete  wall.  There  is  a  possibility  of  a  peak  at  about  18  ms  in  the  experimental  trace 
which  could  be  identified  as  the  secondary  expamion  of  the  explosion.  The  secondary 
expansion  in  the  computation  is  significantly  delayed  due  to  the  gas  bubble 
approximation  and  occurs  at  ~  22  ms. 


Figure  8.  Plot  of  Overpressure  versus  Time  for  Gauge  2  in  Event  3. 


The  pressure  trace  given  in  Figure  9  and  the  peak  pressures  given  in  Table  1  for 
gauge  3  show  reasonable  agreement  between  the  experiment  and  the  calculation.  The 
peak  pressure  for  the  calculation  overestimates  the  experimental  peak  pressure  but 
only  by  ~3  kPa.  The  most  noticeable  difference  between  the  experimental  and 
computational  pressure  traces  is  the  pressure  peak  in  the  experimental  trace  at  about 
23  ms  compart  with  a  dip  in  the  computational  trace  at  that  time.  Given  the  time 
delay  between  the  initial  peak  and  this  peak,  about  15  ms,  this  peak  is  consistent  with 
a  reflection  off  the  angled  end  of  the  adjacent  concrete  wall.  The  absence  of  this  peak 
in  the  calculation  is  probably  due  to  the  approximate  representation  of  this  section  of 
the  concrete  wall  (see  Figure  10).  Since  the  calculation  sees  this  angled  piece  of  wall  in 
a  step-wise  manner,  the  reflection  off  this  piece  of  wall  must  be  lost  due  to  the 
coarseness  of  the  representation. 
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Time  (ms) 


Figure  9.  Plot  of  Overpressure  versus  Time  for  Gauge  3  in  Event  3. 


Figure  10.  Angled  wall  (dashed  line)  and  its  stepwise  representation  in  the  computation  (solid 
line). 


A  further  possibility  for  the  inability  to  model  the  experimental  reflected  pressure 
peak  correctly  could  be  the  orientation  of  the  pressure  gauge.  The  experimental 
pressure  gauge  was  mounted  in  a  baffle  plate  which  was  positioned  such  that  the 
sensing  surface  of  the  gauge  faced  and  was  parallel  to  the  wall.  This  results  in  clean 
traces  when  the  shock  wave  is  travelling  across  the  baffle  plate  (ie,  side  on  or  incident 
pressure)  but  acts  as  a  reflecting  surface  when  the  shock  wave  is  travelling  towards  the 
baffle  plate.  In  this  case,  the  shock  wave  which  has  diffracted  around  the  comer  of  the 
building,  travels  across  the  baffle  plate  resulting  in  a  clean  trace.  However,  the  shock 
wave  which  has  reflected  off  the  retaining  wall  and  then  off  the  wall  of  the  building 
travels  towards  the  baffle  plate,  resulting  in  higher  pressures  due  to  the  interaction 
with  the  baffle  plate.  This  is  represented  schematically  in  figure  11. 
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Figure  11.  Diagram  of  setup  of  pressure  gauges  2, 3  and  4 

The  pressure  trace  given  in  Figure  12  and  the  peak  pressures  given  in  Table  1  for 
gauge  4  show  qualitative  agreement  between  the  experiment  and  the  calculation.  The 
peak  pressure  is  overestimated  but  is  well  within  expectations.  The  reflection  from  the 
concrete  retaining  wall  (at  about  20  ms)  shows  good  time  agreement  but  the 
magnitude  is  substantially  underestimated.  This  is  probably  due  to  the  reasons 
discussed  for  gauge  3.  The  experimental  trace  is  also  complicated  by  a  "ringing",  due 
to  the  baffle  plate  oscillating  at  a  particular  frequency. 


0  t)  20  30  40  50 


Time  (ms) 

Figure  12.  Plot  of  Overpressure  versus  Time  for  Gauge  4  in  Event  3. 
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Figures  13a  -  13h  show  pressure  contour  plots  for  Event  3  at  a  height  of  0.9  m 
above  the  ground,  i.e.,  at  the  height  of  the  midpoint  of  the  explosive  and  the  pressure 
gauges.  Figure  13a  has  an  approximate  outline  of  the  buildings  and  the  wall 
superimposed  on  the  plot  to  give  an  idea  of  their  location.  Unfortunately  the 
structures  are  not  added  onto  the  contour  plots  when  they  are  drawn.  The  contour 
plots  show  the  initial  shock  expanding  out  from  the  explosion  site  (13b),  diffracting 
aroimd  the  comers  of  the  building  (13c,13d)  and  reflecting  off  the  concrete  retaining 
wall  (13e,13f,13g).  The  contour  plots  provide  a  clearer  understanding  of  what  is 
happening  during  the  course  of  the  movement  of  the  shock  front. 
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Event  10. 


The  experimental  results  (solid  line)  and  the  results  of  our  calculations  (dashed  line) 
for  pressure  gauges  1  to  4  in  Event  10  are  given  in  Figures  14  to  17  respectively.  The 
peak  pressures  for  this  event  are  tabulated  in  Table  2.  Note  that  the  calculations  were 
done  with  a  reflection  through  the  axis  of  the  buildings  for  computational  simplicity. 

Table  2.  Peak  experimental  and  calculational  pressures  for  gauges  lto4  in  Event  10. 


Gauge  Number 

Peak  Overpressure  (kPa) 
Experiment 

Peak  Overpressure  (kPa) 
Simulation 

1 

780 

700 

2 

62 

47 

3 

33.6 

23.6 

4 

15.9 

14.9 

It  can  be  seen  from  the  plot  of  the  pressure- time  traces  in  Figure  14  and  the  peak 
pressure  values  given  in  Table  2  that,  for  gauge  1  in  Event  10,  there  is  much  closer 
agreement  than  was  the  case  for  Event  3.  This  is  due  to  the  fact  that  the  pressme 
gauge  was  further  away  from  the  explosion  (6  m  as  compared  to  4  m  in  Event  3)  and 
consequently  the  shock  front  has  longer  (i.e.,  more  computational  cells)  to  resolve 
itself. 


0  20  40  60  80 


Time  (ms) 


Figure  14.  Plot  of  Overpressure  versus  Time  for  Gauge  1  in  Event  10. 
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The  pressure  traces  for  gauge  2  in  Event  10  are  given  in  Figure  15  and  the  peak 
pressures  are  given  in  Table  2.  These  pressure  traces  are  in  good  agreement,  with  the 
only  difference  being  attributable  to  the  secondary  expansion  of  the  explosion 
(at  18  ms  in  the  experimental  trace  and  probably  part  of  the  peak  at  about  30  ms  in  the 
computation). 


Figure  15.  Plot  of  Overpressure  versus  Time  for  Gauge  2  in  Event  10. 


The  pressure  traces  for  gauge  3  in  Event  10  are  given  in  Figure  16  and  the  peak 
pressures  are  given  in  Table  2.  The  traces  show  good  overall  qualitative  agreement, 
however  certain  fine  features  show  a  marked  difference.  There  is  a  small  pressure 
spike  before  the  main  diffracted  shock  reaches  the  gauge.  This  spike,  which  was 
present  in  the  Combustion  Dynamics  2D  cylindrical  calculation[3],  was  identified  as 
the  shock  wave  travelling  through  the  room  in  the  building  and  arriving  at  the  gauge 
before  the  shock  front  which  diffracted  around  the  building.  This  was  not  evident  in 
the  current  calculations  and  is  likely  due  to  the  inclusion  of  the  blast  hole  caused  by 
the  explosive  in  the  Combustion  Dynamics  calculation.  In  the  present  calculation  it 
was  assumed  that  the  structural  response  of  the  building,  i.e.,  the  creation  of  the  blast 
hole,  occurred  on  a  much  longer  timescale  than  the  initial  shock  wave  travelling 
around  the  building. 

The  reflection  of  the  shock  wave  off  the  retaining  wall  arrives  at  approximately  the 
same  time  in  the  experimental  and  computational  traces  (~22  ms)  but  is  much  stronger 
in  the  experimental  case.  This  difference  is  attributed  to  the  orientation  of  the  baffle 
plate  on  the  pressure  gauge  (see  the  discussion  on  gauge  3  in  Event  3).  In  the 
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calculated  pressure  trace  there  is  also  a  large  peak  at  about  50  ms  which  is  not  evident 
in  the  experimental  pressure  trace.  The  contours  in  Figures  18g  and  18h  give  some 
evidence  that  this  is  caused  by  the  second  reflection  of  the  shock  wave  off  the  retaining 
wall  combining  with  the  first  reflection  of  the  secondary  expansion  off  the  retaining 
wall.  In  the  calculation  the  secondary  occurs  later  than  in  the  experiment  due  to  the 
gas  bubble  approximation.  Thus  in  the  experiment  these  multiply  reflected  shock 
waves  show  up  as  separate  shocks  with  considerably  smaller  magnitudes. 


Figure  16.  Plot  of  Overpressure  versus  Time  for  Gauge  3  in  Event  10. 


The  pressure  traces  for  gauge  4  in  Event  10  are  given  in  Figure  17  and  the  peak 
pressures  are  given  in  Table  2.  The  initial  diffracted  peak  in  the  experimental  trace 
shows  good  agreement  with  the  calculated  diffracted  peak.  The  first  reflected  peak 
(about  20  ms)  in  the  calculation  just  precedes  the  experimental  peak  in  time  but  the 
size  of  the  peak  is  in  reasonable  agreement.  The  next  peak  in  the  experimental  trace 
(about  35  ms)  is  missing  from  the  computation  due  to  the  omission  of  the  angled 
section  of  the  concrete  retaining  wall  near  gauge  4  as  a  simplifying  measure  (see 
Figure  2  for  the  experimental  layout).  There  is  good  agreement  with  the  next  reflected 
peak  at  about  45  ms. 
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Figure  17.  Plot  of  Overpressure  versus  Time  for  Gauge  4  in  Event  10. 


The  contour  plots  for  Event  10  are  given  in  Fig;ure  18a  -  18k  and  are  taken  at  a 
height  of  0.7  m.  The  development  of  the  shock  wave  expanding  out  from  the 
explosion  site  (18b,  18c),  reflecting  off  the  concrete  wall  (18d,  18e,  18f),  diffracting 
around  the  comer  of  the  building  (18c-18g)  and  inside  the  room  in  the  building  (18d- 
18i)  can  all  be  followed  using  these  contour  plots.  The  outstanding  difference  between 
event  3  and  event  10  is  the  inclusion  of  the  internal  room  in  event  10  which  is  clearly 
shown  in  the  contour  plots. 
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Figure  18i.  563ms  contours 


Figure  18 j.  67.1ms  contours 


5.  Conclusions 


This  paper  represents  the  first  successful  modelling  of  the  Australian  Challenge  in  full 
3-D  geometry.  The  results  in  this  paper  have  shown  that  the  program  has  the  ability 
to  successfully  calculate  the  blast  overpressure  from  an  explosion  in  a  structurally 
complex  scenario.  This  validation  of  the  code  will  allow  it  to  be  used  with  confidence 
to  model  complex  operational  scenarios  requiring  the  rapid  entry  of  forces  into 
buildings.  It  will  reduce  dramatically  the  need  for  expensive  and  time-consuming 
trials  in  this  area  and  allow  hypothetical  situations  to  be  explored.  It  will  also  allow 
the  safety  of  current  operational  procedures  to  be  evaluated,  as  well  as  any  proposed 
changes  to  these  procedures. 

There  are  certain  limitations  on  the  capabilities  of  the  code  for  doing  these  types  of 
calculations  given  the  size  of  the  currently  available  computers.  First,  a  computer 
memory  size  limitation  means  that  a  grid  of  approximately  1  million  cells  is  the 
maximum  that  can  be  considered.  This  in  turn  directly  limits  the  fineness  of  the  grid 
that  can  be  considered,  for  example,  for  a  grid  of  10  m  x  10  m  x  10  m  the  minimum  cell 
size  available  is  10  cm  x  10  cm  x  10  cm.  A  less  significant  limitation  is  the  duration  of 
the  calculation  which  is  currently  restricted  to  about  3000  timesteps.  This  allows  the 
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calculation  to  be  completed  within  about  1  week.  This  limitation  precludes  the  use  of 
the  code  in  a  real-time  operational  situation  -  the  results  will  simply  take  too  long  to  be 
of  operational  use. 

A  further  limitation  is  that  the  program  employs  a  Cartesian  grid  and 
consequently  can  only  approximately  represent  non-Cartesian  structures  in  the 
calculation.  An  example  of  this  in  the  current  calculations  is  the  bent  part  of  the  1.8  m 
high  wall  near  the  buildings.  Since  this  part  of  the  wall  is  at  a  non-90°  angle  it  is 
represented  by  a  stepwise  construction  rather  than  a  smooth  piece  of  wall.  Similarly 
circular  structures  would  not  be  well  represented  by  the  Cartesian  gnd  because  of  the 
limitations  imposed  by  the  resolution  of  the  calculation. 

These  limitations  on  the  calculations  do  not  present  major  restrictions  on  the  use  of 
the  program  but  should  be  kept  in  mind  when  deciding  whether  this  program  will 
give  the  accuracy  required  in  a  desired  situation.  This  work  has  shown  that  this 
program  can  successfully  calculate  the  blast  overpressure  from  an  explosion  in  a  large 
complex  scenario  of  interest  to  the  Australian  Defence  Forces.  Consequently  any 
giimilar  problem,  for  example,  the  rapid  entry  of  the  Special  Forces  into  buildings,  will 
be  able  to  be  modelled  with  confidence  and  without  resort  to  expensive  field  trials. 
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Appendix  1: 

Inclusion  of  an  Arbitrary  Geometry  in  the  Program 

The  program  is  written  in  FORTRAN  and  so  all  of  the  programming  examples  iised  in 
this  appendix  use  a  FORTRAN  style  to  illustrate  the  point. 

The  program  has  a  3-dimensional  Cartesian  grid  as  a  basis  for  its  computations. 
That  is,  when  the  program  is  integrating  in  each  of  the  three  directions  it  will  loop  over 
the  number  of  points  in  each  direction,  e.g., 

do  100  k=l,ncz 
do  110  j=l,ncy 
do  120  i=l,ncx 


120  continue 
110  continue 
100  continue 

In  this  case  all  the  cells  are  being  calculated  since  they  are  all  included  in  the  do- 
loop.  To  illustrate  the  way  in  which  arbitrary  geometries  are  included  in  three 
dimensions  we  will  examine  the  2-dimensional  case.  We  then  have  the  do-loops 

do  100  j=l,ncy 
do  110  i=l,ncx 


110  continue 
100  continue 

and  this  can  be  represented  graphically  as 
NCYl - 


II _ I 

1  NCX 
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If  we  now  want  to  exclude  part  of  the  area,  e.g.,  between  x=ncx/3  and  x=2ncx/3 
then  we  need  to  modify  the  do-loops  to  be 

do  100  j=l,ncy 
do  110  i=l,ncx/3 
110  continue 

do  120  i=2ncx/3,ncx 


120  continue 
100  continue 


1  NCX/3  2NCX/3  NCX 


This  approach  can  be  generalised  to  an  arbitrary  number  of  excluded  areas 

do  100  j=l,ncy 
do  110  iarea=l,inumarea 
ixmin=ileft(iarea) 
ixmax=iright(iarea) 
do  120  i=ixmin,ixmax 


120  continue 
110  continue 
100  continue 

where  inumarea  is  the  number  of  areas  to  be  integrated,  ixmin  is  the  left  limit  for  that 
area  and  ixmax  is  the  right  limit  for  that  area.  If  further  areas  were  to  be  excluded  but 
which  only  extended  part  of  the  way  along  the  y  axis,  e.g.. 
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then  this  can  be  done  in  the  following  way 

do  100  j=l,ncy 
ixnumber=ixstep(i) 
do  110  istep=l,ixnumber 
ixmin=ileft(j4xnumber) 
banax=iright(j,ixnumber) 
do  120  i=ixmin,ixmax 


120  continue 
110  continue 
100  continue 

where  the  array  ixstep  is  the  number  of  integrations  for  the  value  of  j  and  the  arrays 
ileft  and  iright  are  the  left  and  right  limits  for  each  integration. 

The  next  step  is  how  to  quickly  and  easily  enter  the  limits  of  the  integration 
represented  in  the  arrays  ileft  and  iright.  The  method  chosen  in  this  approach  can  be 
easily  shown  graphically  below  where  the  final  geometry  is  simply  built  up  by  a  series 
of  overlaying  levels  which  contain  a  rectangle,  either  solid  to  make  an  obstruction  or 
empty  to  open  up  a  space. 
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NCY 


1 

1  NCX/3  2NCX/3  NCX 


Addition  of  Solidl 


Addition  of  Solid2 


Addition  of  Openl 


Addition  of  Open2 


To  take  a  more  specific  example,  if  we  let  ncx  =  ncy  =  30  then  we  can  specify  the 
rectangles  by  their  beginning  and  end  points  along  each  axis.  For  the  solid  rectangles 
we  have  the  coordinates 


xmin 

xmax 

ymin 

ymax 

solidl 

10 

20 

1 

30 

solid2 

1 

10 

22 

30 

and  for  the  open  rectangles  we  have  the  coordinates 
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xmin  xmax  ymin  ymax 

openl  10  20  14  17 

open2  15  20  2  8 

where  solidl  has  the  comers  (proceeding  clockwise  from  the  lower  left) 
(10,1),(10,30),(20^0),(20,1). 

The  question  remains  of  how  do  we  translate  these  numbers  to  obtain  the  correct 
limits  for  the  integrations.  Again  this  can  best  be  shown  by  example.  We  shall 
examine  several  specific  points  along  the  y  axis,  namely,  y=10,  y=28,  y=15  and  y=5 
which  will  show  the  addition  of  each  rectangle  (these  hnes  are  shown  in  the  figure 
below) 


1  NCX/3  2NCX/3  NCX 


For  y=10  the  integration  limits  after  each  addition  are 


left 

right 

initial 

1 

30 

additions 

solidl 

1 

10 

20 

30 

solid2 

1 

10 

20 

30 

openl 

1 

10 

20 

30 

open  2 

1 

10 

20 

30 

final 

1 

10 

20 

30 
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For  y=28  the  integration  limits  after  each  addition  are 


left 

right 

initial 

1 

30 

additions 

solidl 

1 

10 

20 

30 

solid2 

20 

30 

openl 

20 

30 

open  2 

20 

30 

final 

20 

30 

For  y=15  the  integration  limits  after  each  addition 

left 

right 

initial 

1 

30 

additions 

solidl 

1 

10 

20 

30 

solid2 

1 

10 

20 

30 

openl 

1 

30 

open  2 

1 

30 

final 

1 

30 

For  y=10  the  integration  limits  after  each  addition 

left 

right 

initial 

1 

30 

additions 

solidl 

1 

10 

20 

30 

solid2 

1 

10 

20 

30 

openl 

1 

10 

20 

30 

open  2 

1 

10 

15 

30 

final 

1 

10 

15 

30 
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That  is  the  limits  for  a  particular  y  value  change  whenever  the  rectangle  includes 
that  value.  In  this  way  we  can  build  up  quite  complex  shapes  by  adding  either  solid 
or  open  rectangles  and  consequently  adjusting  the  limits  of  integration  for  each 
addition. 

These  concepts  and  procedures  are  quite  easily  extended  to  three  dimensions  by 
adding  a  further  do-loop,  e.g., 

do  100  k=l,nc2 
do  110  j=l,ncy 
ixnumber=ixstep(j,k) 
do  120  istep=l,ixnumber 
ixmin=ileft(j,k,ixnumber) 
ixmax=iright(j,k,ixnumber) 
do  130  i=ixmin,ixmax 


130  continue 
120  continue 
110  continue 
100  continue 
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